﻿/*
 * SphericConformalTransform
 * 
 * initial implementation : 6.7.2011
 *
 * Copyright (c) Adam Sida, asida@centrum.cz.
 * All rights reserved.
 * 
 * for license see Petra.Gis.License.txt
 */
using System;
using System.Collections.Generic;
using System.Runtime.Serialization;
using Petra.Gis.Geometry.Support;


namespace Petra.Gis.ReferenceSystemServices.Transformations
{

  /// <summary>
  /// This is pseudo-datum transform especially used in european coordinate systems
  /// (Swiss, Hungary, Czech, Slovak) where standard cartographic projection (conic, cylindrical)
  /// is at general (oblique) position on spheroid. This projection is one more step 
  /// that must be done in transform between geographic space (lambda, phi) and specific projection (x,y).
  /// Projection uses gauss sphere with radius given by specific standard parallel.
  /// Sphere position relative to Spheriod is also given by oblique standard parallel,....
  /// This projection supports also pole reprojection (cartographic pole). 
  /// New pole position is given by its position in gauss spheric space (lambda, phi).
  /// Transform can threfore produce not only standard (lambda, phi) coordinates, but also
  /// non standard (rotated) spheric coordinates. Coordinate names are defined to U (longitude), V (lattitude).
  /// Projection has different equations for forward and reverse direction.
  /// Forward direction is defined as gauss sphere to ellipsoid.
  /// </summary>
  [DataContract]
  public class SphericConformalTransform : InversibleParametricTransform
  {

    //basic transform constants
    private double Phi0;  // primary (standard) parallel (undistorted = on spheriod)
    private double alpha; // projection constant
    private double k;     // projection constant k
    private double ka;    // inverse to k
    private double e;     // eccentricity
    private double eHalf; // e/2

    //private const int MAXIMUM_ITERATIONS = 15;
    private const double ITERATION_TOLERANCE = 1E-11;

    //equation variables
    private double tmpLat;
    private double sinPhi;
    private double cosPhi;
    private double esf;
    private double fi1;
    private double kaU;
    private double dV;
    //private double U0; //primary (standard) parallel (undistorted = on gauss sphere)
    
    //oblique constants
    private bool isOblique;
    private double Uc; //cartographic pole on gauss sphere (phi)
    private double Vc; //cartographic pole on gauss sphere (lambda)
    private double sinUc;
    private double cosUc;
    

    //private double sphereX;
    //private double sphereY;
    //private double dfy;

    AngularUnit srcUnit;
    AngularUnit trgtUnit;




    #region constructors

    /// <summary>
    /// xaml constr.
    /// </summary>
    public SphericConformalTransform() 
      :base()
    {
      Metadata = new Metadata("SphericConformalTransform", "none", -1, "", "", "provides transformation between Geographics Coordinate systems and gauss sphere");
    }





		public SphericConformalTransform(GeographicCoordinateSystem targetCS,
                                     double obliqueStandardParallel,
                                     AngularUnit angularParamUnit,
                                     bool isInverse)
			:this(targetCS, angularParamUnit, isInverse,
            new TransformParameters(new XParam("oblique_standard_parallel", obliqueStandardParallel)) )
    {  }





    protected SphericConformalTransform(GeographicCoordinateSystem targetCS,
                                       AngularUnit angularParamUnit,
                                       bool createInverse,
                                       TransformParameters parameters)
      : base(GeographicCoordinateSystem.CreateGaussSphereCS(targetCS, parameters["oblique_standard_parallel"] * angularParamUnit.RadiansPerUnit),
             targetCS, null, angularParamUnit, createInverse, parameters)
    {
      if (parameters.Count < 1)
        throw new ArgumentException("SphericConformalTransform requires at least 1 parameter");

      if(targetCS.Unit == null)
        throw new ArgumentException("coordinate system must have unit defined");

      Metadata = new Metadata("SphericConformalTransform", "none", -1, "", "", "provides transformation between Geographics Coordinate systems and gauss sphere");
    }


    #endregion constructors





    protected override void updateTransform()
    {
      //Check for missing parameters
      //if (!Parameters.ContainsKey("gauss_sphere_standard_parallel_1"))
      //  throw new ArgumentException("Missing projection parameter 'gauss_sphere_standard_parallel_1'");

      srcUnit = (AngularUnit)Source.Unit;
      trgtUnit = (AngularUnit)Target.Unit;

      Phi0 = Parameters["oblique_standard_parallel"] * base.AngularParameterUnit.RadiansPerUnit;

      //optional parameters, cartgraphic pole
      //given as geog. position on gauss sphere !!!
      if (Parameters.ContainsKey("oblique_pole_latitude"))
        Uc = Parameters["oblique_pole_latitude"] * base.AngularParameterUnit.RadiansPerUnit;

      if (Parameters.ContainsKey("oblique_pole_longitude"))
        Vc = Parameters["oblique_pole_longitude"] * base.AngularParameterUnit.RadiansPerUnit;

      //eccentriticy
      e = ((GeographicCoordinateSystem)Target).Datum.Ellipsoid.Eccentricity;
      eHalf = e / 2.0;

      //square of semi-major axis
      double e2 = ((GeographicCoordinateSystem)Target).Datum.Ellipsoid.e2;

      //semi-major axis
      //double a = Projection.GeographicCoordinateSystem.Datum.Ellipsoid.SemiMajorAxis
      //  *Projection.GeographicCoordinateSystem.Datum.Ellipsoid.AxisUnit.MetersPerUnit;

      //const alpha
      alpha = Math.Sqrt(1 + (e2 * Math.Pow(Math.Cos(Phi0), 4) / (1 - e2)));

      //primary (standard) parallel (undistorted = on gauss sphere)
      double U0 = Math.Asin(Math.Sin(Phi0) / alpha);

      //const k
      k = Math.Pow(Math.Tan(Phi0 / 2 + ProjMath.S45), alpha)
        * Math.Pow((1 - e * Math.Sin(Phi0)) / (1 + e * Math.Sin(Phi0)), (alpha * e / 2.0))
        / Math.Tan(U0 / 2 + ProjMath.S45);

      ka = Math.Pow(1.0 / k, -1.0 / alpha);

      ////first vertical curvature radius of spheroid
      //double N = a / (Math.Sqrt(1 - e2 * Math.Pow(Math.Sin(Phi0), 2)));
      //double M = a * (1 - e2) / Math.Pow(1 - e2 * Math.Pow(Math.Sin(Phi0), 2), 2 / 3);

      ////gauss sphere radius
      //double R = Math.Sqrt(M * N);


      if (Uc != 90 || Vc != 0)
      {
        sinUc = Math.Sin(Uc);
        cosUc = Math.Cos(Uc);
        isOblique = true;
      }
      else
        isOblique = false;

    }





    protected override void transform(ICoordinate coord)
    {

      if (IsInverse)
      {
        //from datum to aposphere/gauss sphere...
        coord.lon *= trgtUnit.RadiansPerUnit;
        coord.lat *= trgtUnit.RadiansPerUnit;

        sinPhi = Math.Sin(coord.lat);

        //equation verified and optimized by Kristyna Ivankova
        coord.lat = 2 * Math.Atan(1 / k * Math.Pow(Math.Sqrt((1 + sinPhi) / (1 - sinPhi))
          * Math.Pow((1 - e * sinPhi) / (1 + e * sinPhi), eHalf), alpha)) - (ProjMath.HALF_PI);

        coord.lon = coord.lon * alpha;

        if (isOblique)
        {
          // spheric triangulation to cartogr. pole on sphere...
          // U,V -> karto
          //lat
          tmpLat = Math.Asin(sinUc * Math.Sin(coord.lat) + cosUc * Math.Cos(coord.lat) * Math.Cos(Vc - coord.lon));
          //lon
          coord.lon = Math.Asin(Math.Sin(Vc - coord.lon) * Math.Cos(coord.lat) / Math.Cos(tmpLat));

          coord.lat = tmpLat;
        }

        //return coordinate uses same units as spheroid
        coord.lon /= srcUnit.RadiansPerUnit;
        coord.lat /= srcUnit.RadiansPerUnit;

      }
      else 
      {
        //from aposphere/gauss sphere to datum
        coord.lon *= srcUnit.RadiansPerUnit;
        coord.lat *= srcUnit.RadiansPerUnit;

        if (isOblique)
        {
          //spheric triangulation karto -> U,V
          tmpLat = Math.Asin(Math.Sin(coord.lat) * sinUc - Math.Cos(coord.lat) * cosUc * Math.Cos(coord.lon));//U
          coord.lon = Vc - Math.Asin(Math.Cos(coord.lat) * Math.Sin(coord.lon) / Math.Cos(tmpLat)); //V
          coord.lat = tmpLat;
        }

        //throw new NotImplementedException();

        //double ro = ProjMath.Hypot(coord., y);
        //double eps = Math.atan2(-x, -y);
        //double d = eps / n;
        //double s = 2 * (Math.Atan(Math.Pow(ro0 / ro, 1 / n) * tanS2) - s45);

        cosPhi = Math.Cos(coord.lat);
        //double u = Math.Asin((cosAzim * Math.Sin(coord.lat)) - (sinAzim * cosPhi * Math.Cos(coord.lon)));
        kaU = ka * Math.Pow(Math.Tan((coord.lat / 2.0) + ProjMath.S45), 1.0 / alpha);
        dV = Math.Asin((cosPhi * Math.Sin(coord.lon)) / Math.Cos(coord.lat));

        //coord.lon = -dV / alpha;
        coord.lon = dV / alpha; //todo:might need update !!!

        //resolve lat...
        tmpLat = 0;
        fi1 = coord.lat;

        // iteration calculation
        for (int i = ProjMath.MAX_ITERATIONS;;)
        {
          fi1 = tmpLat;
            esf = e * Math.Sin(fi1);
            tmpLat = 2.0 * (Math.Atan(kaU * Math.Pow((1.0 + esf) / (1.0 - esf), eHalf)) - ProjMath.S45);
            if (Math.Abs(fi1 - tmpLat) <= ITERATION_TOLERANCE)
                break;

            if (--i < 0) 
                throw new ApplicationException();

        }

        coord.lat = tmpLat;

        coord.lon /= trgtUnit.RadiansPerUnit;
        coord.lat /= trgtUnit.RadiansPerUnit;
        
      }
    }





    protected override bool parametersValid()
    {
      return Parameters.ContainsKey("oblique_standard_parallel");
    }





    /// <summary>
    /// must override here, spheric conformal doesnt need parameter distance units (only angular)
    /// </summary>
    /// <returns></returns>
    protected override bool graphCompleted()
    {
      return !(this.AngularParameterUnit == null || !parametersValid());
    }





    public override SRSTransform Inverse()
    {
      return new SphericConformalTransform(
                      (GeographicCoordinateSystem)Target,
                      AngularParameterUnit,
                      !IsInverse,
                      Parameters);
    }





    public override Dimension TransformDimension
    {
      get { return Dimension.dim2D; }
    }





    /// <summary>
    /// creates new instance of Spheric conformal transform from geogr. coordinate system
    /// </summary>
    /// <param name="gcs">source georg. coordinate system</param>
    /// <param name="angParamUnit">transformation parameter angular unit (distance is not relevant for SphericConformal)</param>
    /// <param name="parameters">transform parameters</param>
    /// <returns>new instance of SphericConformalTransform or null or throws exception when invalid input parameters</returns>
    public static SphericConformalTransform Create(GeographicCoordinateSystem gcs, AngularUnit angParamUnit, bool createInverse, TransformParameters parameters)
    {
      if (parameters["oblique_pole_latitude"] * angParamUnit.RadiansPerUnit > ProjMath.HALF_PI)
        throw new ArgumentOutOfRangeException("latitude must be in 0 - 90 deg range");

      return new SphericConformalTransform(gcs, angParamUnit, createInverse, parameters);
    }


  }

}
